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Abstract 


A  complex  variable  boundary  element  method  is  developed  for  potential  flow 
problems  by  applying  Cauchy's  integral  theorem  to  the  complex  velocity.  The 
resulting  integral  equation  is  a  function  of  the  normal  and  tangential  velocity 
components  on  the  boundary.  A  new  form  of  the  full  nonlinear  dynamic  free  surface 
boundary  condition  is  used  to  describe  the  evolution  of  tangential  velocities.  This 
alternate  method  solves  for  flows  with  field  singularities  more  easily  than  the 
conventional  method,  which  uses  the  complex  velocity  potential.  Also,  the  velocity 
rield  is  given  directly  without  the  need  for  numerical  differentiation.  Under  the  new 
formulation,  the  dynamic  free  surface  boundary  condition  does,  however,  become 
more  complicated.  As  a  result,  while  the  new  form  of  the  boundary  element  method 
has  definite  advantages  for  fixed  boundaries,  its  usefulness  for  free  surface  problems 
is  mixed. 
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1  Introduction 


Since  Longuet-Higgins  and  Cokelet  [1]  introduced  a  semi-Lagrangian  method  for  the  full 
nonlinear  initial-boundary  value  problem  of  breaking  wave  simulation,  a  number  of  different  nu¬ 
merical  models  for  fully  nonlinear  free  surface  problems  have  been  developed,  e.g.,  Vinje  and 
Brevig  [2];  Baker,  Meiron  and  Orszag  [3];  Dold  and  Peregrine  [4];  etc.  Longuet  Higgins  and 
Cokelet  [1]  used  Green’s  theorem  to  develop  an  integral  equation  for  the  normal  velocity,  §£, 
when  the  velocity  potential,  <p ,  is  given  as  an  initial  condition  on  the  free  surface.  Vinje  and 
Brevig  [2]  applied  Cauchy’s  integral  theorem  to  the  comp'ex  potential  (0  =  <j>  +  tip)  to  form  an 
integral  equation  for  ip.  Baker,  Meiron,  and  Orszag  [3]  developed  a  method  wherein  the  velocity 
field  is  represented  by  distributed  vortices  (or  dipoles)  along  the  free  surface.  The  time  evolution 
of  the  vortex  sheet  strength  is  found  by  solving  an  integral  equation.  In  all  these  cases,  numerical 
differentiation  or  integration  is  needed  to  obtain  the  velocity  on  the  free  surface. 

Martensen  [5]  and  Lewis  and  Ryan  [6]  employed  boundary  element  methods  directly  to  deter¬ 
mine  the  velocity  on  boundaries.  Their  work  was  not  extended  to  free  surface  problems.  Recently, 
Dold  and  Peregrine  [4]  used  a  boundary  element  method  to  directly  compute  the  complex  veloc¬ 
ity  (k)  based  on  Cauchy’s  integral  theorem.  Their  form  of  the  dynamic  free  surface  boundary 
condition  is  based  on  the  velocity  potential,  so  recalculation  of  the  rate  of  change  of  the  velocity 
potential  following  a  material  particle  was  necessary. 

In  real  flows,  viscosity  generates  vorticity  at  the  surface  of  solid  bodies.  This  vorticity  is  shed, 
resulting  in  an  unsteady  wake  flow.  Vortex  methods  (e.g.,  [7])  assume  that  vorticity  is  shed  as 
discrete  vortices  or  as  a  vortex  sheet  in  an  otherwise  irrotational  flow.  In  the  two-dimensional 
problem,  the  velocity  potential  of  the  shed  vortices  is  multivalued  (i.e.,  the  domain  has  branch 
cuts).  Special  care  and  treatment  is  therefore  needed  to  incorporate  these  vortices  into  velocity 
potential  boundary  element  methods. 
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An  alternative  two-dimensional  boundary  element  model  for  complex  velocity  is  derived  here, 
together  with  a  modified  dynamic  boundary  condition  on  the  free  surface  in  terms  of  the  velocity. 
The  advantages  and  disadvantages  are  discussed,  and  some  simple  computational  examples  are 
given. 


2  Formulation  using  Cauchy’s  Theorem 


Usually  Cauchy’s  theorem  is  applied  to  the  complex  potential  to  derive  a  complex  variable 
boundary  element  method  for  irrotational  flow  problems.  Here,  we  use  Cauchy’s  theorem  to 
derive  an  integral  equation  in  terms  of  the  complex  velocity.  Let  the  complex  velocity  of  the 
problem  be  x,  which  is  analytic  inside  the  domain  and  is  given  as 

k  —  u  —  iv  ,  (1) 

where  u  and  v  are  the  velocity  components  in  the  x  and  y  directions,  respectively. 

Figure  1  shows  the  problem  domain  and  bounding  surfaces.  Cauchy’s  theorem  gives 

I  _?(-)  dz  _  ,QK(£t)  >  (2) 

Jan  2  —  U 

where  a  is  0  or  2jt  when  the  location  of  the  kernel  singularity,  £*,  is  outside  or  inside  the  boundary, 
respectively.  If  the  kernel  singularity  is  on  the  boundary  (C*  £  <90),  a  is  equal  to  the  included 
angle,  and  the  integral  is  treated  as  principal- valued.  We  take  C*  to  approach  the  boundary 
from  the  outside  of  the  domain  so  that  a  is  zero.  The  algebraic  system  is  formed  from  (2)  by 
discretizing  the  integral  and  letting  the  kernel  singularity  approach  each  of  the  N  nodal  points 
in  turn,  i.e.  Ok  — ' ►  A  special  limiting  process  is  needed  to  evaluate  the  integral  near  Ct- 
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The  boundary  contour,  3fl,  is  decomposed  into  dSln  and  dCls,  where  the  normal  velocity  is 
given  on  and  the  tangential  velocity  is  given  on  dCl, .  For  a  free  surface  time  marching 
problem,  <90  n  would  represent  a  solid  surface  and  dQ,  a  free  surface.  Hence,  the  boundary 
conditions  can  be  categorized  by  the  following  given  conditions: 


Re{Keie}  =  V, 

on  3@r  , 

(3) 

Im{Ke,e}  =  Vn 

on  dfei  . 

(4) 

Here  9  is  the  tangent  angle  with  respect  to  the  positive  x  axis  as  shown  in  Figure  1. 


The  boundary  contour  is  represented  as  piecewise-linear  panels  and  the  complex  velocity  as 
piece  wise- linear  functions  in  a  panel,  zj  and  kj  denote  the  location  and  complex  velocity  at  the 
jth  node  point.  Our  numerical  model  assumes  that  k  on  the  jth  panel  is  of  the  form 


The  total  complex  velocity  can  be  written  as 

K  =  K[  +  Kv  +  Kd  ,  (9) 

where  subscripts  I,  v,  and  d  denote  the  components  due  to  the  incoming  flow,  the  free  vortices, 
and  the  remainder  (i.e.,  body  and  free  surface  singularities),  respectively.  Usually  in  a  problem 
ki  and  Ky  will  be  given,  and  k&  will  be  unknown. 
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The  complex  velocity  due  to  vortices  is  given  by 

N. 


(10) 


where  y j  is  the  strength  of  a  vortex  located  at  Zj ,  and  Nv  is  the  total  number  of  vortices.  Note 
that  the  complex  velocity  of  a  vortex  is  single-valued,  but  its  complex  potential  is  not.  The 
complex  velocity  potential  of  a  vortex  or  source  is  given  as  a  logarithmic  function,  the  imaginary 
part  of  which  is  multi-valued.  Hence,  in  the  conventional  potential  method  we  have  difficulty  in 
uniquely  satisfying  the  boundary  conditions  due  to  branch  cuts.  This  is  also  true  for  the  integral 
equation  formulation  using  Green’s  theorem. 


We  apply  Cauchy’s  theorem  to  nd  instead  of  k,  since  Kd  is  expected  to  go  to  zero  far  from  the 
region  of  disturbance.  We  could  instead  apply  the  residue  theorem  to  k d  +  kv,  but  the  present 
method  is  computationally  equivalent.  In  either  case,  the  contour  integration  contribution  in 
the  far  field  will  be  eliminated.  The  disturbance  from  a  solid  body  dies  out  as  0(l/z2),  but  the 
disturbance  at  a  free  surface  does  not  decay,  due  to  travelling  waves.  To  make  the  computational 
domain  finite,  either  spatial  periodicity  is  assumed  or  an  initial-value  problem  is  studied  for 
small  time  such  that  the  free  surface  disturbances  are  confined  to  a  small  region.  Alternatively, 
damping-layers  (Baker  et  al.  [8])  can  be  used  to  artificially  dampen  the  far  field  effects  on  the 
free  surface. 

The  total  complex  velocity  on  the  contour  can  be  rewritten  in  terms  of  the  total  normal  and 
tangential  velocity  as 

k  eie  =  V,  +  iV„  .  (II) 

Rewriting  equation  (6)  in  terms  of  Kd  gives  the  following  algebraic  equation  in  terms  of  Vd,  and 

vdn- 

N 

+  *  V^)jRjk  =0  for  k  =  1, . . . ,  N  ,  (12) 

*  i=i 

where  Rjt  —  Tjte~*9>,  and  Vd,  and  Vdn  are  the  tangential  and  normal  disturbance  velocity 
components,  respectively.  Either  Vd,  or  Vd„  is  specified  on  the  boundary  contour  by  the  boundary 
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conditions. 

The  above  set  of  complex  algebraic  equations  gives  twice  as  many  real  equations  as  unknowns. 
It  is  sufficient  to  solve  half  of  the  set  of  equations,  say,  just  the  real  part  of  (12).  Instead,  we 
elect  to  solve  all  of  the  equations  in  a  least-squares  sense  with  the  increased  accuracy  shown  in 
Schultz  and  Hong  [9]. 


3  Free  Surface  Boundary  Conditions 


The  kinematic  and  dynamic  boundary  conditions  of  a  free  surface  for  an  inviscid  flow  are 
given  as 

i 


D<{> 

Dt 


Dz 
Dt 
1 


=  K 


=  -*»  +  §  |V*  |2 


(13) 

(14) 


where  ’  denotes  the  complex  conjugate,  jjj  represents  a  material  derivative,  and  g  is  the  accel¬ 
eration  due  to  gravity  (acting  here  in  the  —y  direction).  The  kinematic  condition  assures  that 
material  particles  remain  on  a  free  surface.  The  dynamic  boundary  condition  is  the  unsteady 
Bernoulli  equation,  requiring  that  pressure  be  constant  along  a  free  surface. 


Since  our  boundary  integral  formulation  is  written  in  terms  of  the  velocity,  the  dynamic 
boundary  condition  should  be  as  well.  The  Euler  equation  can  be  written  as 


DV 


Vp 


n,  - - 93 

Dt  p 


Since  the  pressure  gradient  along  the  free  surface  is  zero,  multiplication  by  the  unit  tangent 
vector  t,  gives 


DV.  v  Dtn 
Dt  "  Dt 


The  first  term  on  the  right-hand  side  results  from  the  local  coordinate  system  being  noninertial. 


Rewriting  this  equation  in  terms  of  9  gives 


DV,  v  DO 

Dt  n  Dt 


gsmO  . 


(15) 


Thus,  the  rate  of  change  in  the  tangential  velocity  following  a  material  particle  on  the  free  surface 
is  given  by  (15)  as  a  function  of  the  normal  velocity,  the  tangent  angle  9,  and  the  rate  of  change 
of  9  following  a  material  particle  (^f).  Since  there  is  a  time  derivative  term  on  the  right-hand 
side,  equation  (15)  is  not  in  a  standard  form  for  applying  a  time  marching  scheme. 

x 

The  material  derivative  of  the  tangent  angle,  £*,  can  be  calculated  from  the  normal  and 
tangential  velocity  distribution  along  a  free  surface  by  (Hong,  TVyggvason,  and  Schultz  [10]) 


Di  =  _dVn 

Dt  ds  +  R 


(16) 


Here,  R  is  the  radius  of  curvature  of  the  boundary  contour  and  is  defined  as  positive  when  the 
center  of  curvature  is  located  on  the  left-hand  side  as  we  follow  the  contour  s. 


Equation  (16)  has  the  disadvantage  of  requiring  two  spatial  derivatives  of  the  contour  to  find 
the  radius  of  curvature.  For  the  highly  deformed  free  surface  contour,  the  numerical  error  in 
calculating  the  curvature  can  be  very  large,  unless  a  high-order  approximation  is  used.  Our  free 
surface  computations  using  (15)  and  (16)  were  satisfactory  only  when  the  free  surface  curvatures 
were  small.  Here,  we  avoid  the  difficulty  of  performing  two  spatial  derivatives  to  find  R,  by  using 
another  approximation  scheme  described  below. 

Since  the  present  method  uses  the  tangential  and  normal  velocity  components,  the  definitions 


of  the  tangent  angle,  9,  and  the  time  rate  of  change  of  the  tangent  angle,  ,  are  crucial  to  the 
overall  accuracy.  We  approximate  the  tangent  angle  at  the  jth  node  9j  as 


-f  e'l‘  Dj-i 
+e'#> 


(17) 


where 
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and  Dj  is  the  length  of  the  jth  panel.  Hence,  the  tangent  angle  at  a  node  is  given  by  a  weighted 
average  of  the  mean  tangent  angle  of  the  adjacent  panels. 


A  similar  procedure  is  taken  to  approximate  the  time  rate  of  change  of  the  tangent  angle. 
at  node  j  is  approximated  by  the  weighted  average  of  8j-i  and  Oj  as 


where 


DOj  _  0)  - 1 D j  +  Oj  Dj  - 1 

Dt  Dj  ■+•  -Dj-i 


0 i  /w 

l  -  h  J 


(19) 


is  the  time  derivative  of  Oj  given  by  (18).  We  find  (19)  to  be  much  more  accurate  than  (16)  when 
piecewise  linear -panels  are  used. 

4  Numerical  Calculation  and  Examples 


The  algorithm  is  tested  on  three  examples  to  demonstrate  its  usefulness  and  accuracy.  Fol¬ 
lowing  Longuet-Higgins  and  Cokelet  [1],  an  initial-boundary  value  problem  is  defined  for  each 
case  such  that  V,  is  given  on  a  free  surface  as  an  initial  condition  and  V„  is  given  on  the  solid 
body  surface. 


We  use  a  fourth-order  Runge-Kutta-Gill  scheme  to  evolve  the  free  surface  variables.  The 
overdetermined  matrix  problem  is  solved  using  an  iterative  conjugate  gradient  method  [11],  This 
solver  is  very  efficient  when  a  good  initial  guess  is  available,  aw  in  time  marching  problems. 

The  first  example  is  a  numerical  simulation  of  a  propagating,  nonlinear  gravity  wave.  The 
results  of  the  present  method  are  compared  with  the  potential  method.  The  second  example  is 
a  vortex  pair  and  wall  interaction  problem.  The  last  example  studies  a  vortex  pair  interacting 
with  a  free  surface.  For  the  second  and  last  example,  no  comparisons  are  made  with  the  poten¬ 
tial  method  due  to  the  expected  difficulties  of  using  the  potential  method  for  the  field  vortices 
problem,  as  mentioned  in  section  2. 


4.1  Nonlinear  Wave  Simulation 

A  deep  water,  periodic,  nonlinear  gravity  wave  problem  is  solved  using  the  present  method 
and  the  results  are  compared  with  those  of  the  complex  potential  boundary  element  method 
(Schultz  [12]). 

The  initial  surface  elevation,  complex  potential,  and  corresponding  tangential  velocity  are 
given  in  Table  1.  Initial  nodes  are  spaced  uniformly  in  the  x  direction.  Since  this  problem  is 
periodic  in  x,  the  computational  domain  is  reduced  to  a  finite  region  in  the  x  direction  using 
conformal  mapping  [12] 


c  =  exp(—2iri  c/A),  (20) 

where  A  is  the  wavelength. 

Error  distributions  in  the  x  and  y  velocity  components  at  the  initial  time  are  compared  for  the 
two  methods  in  Figures  2a  and  2b.  A  central  difference  scheme  is  used  to  calculate  the  velocity 
in  the  potential  formulation.  The  velocity  formulation  is  seen  to  give  slightly  better  results,  but 
the  accuracy  of  both  methods  is  consistent  with  the  precision  of  the  calculation. 

In  Figures  3a,  b,  c,  and  d,  the  time  variation  of  the  maximum  surface  elevation  (ym4r),  the 
minimum  surface  elevation  (ymjn),  the  mass  conservation  (Em),  and  the  total  energy  (E(t))  is 
shown  for  both  methods.  Here  Em  and  E(t)  are  defined  by 

Em  =  J  n(x,t)dx/i^J  rj2(x,0)d;rj  (21) 

m  =  \j  *7 2(x,t)  dx  +  £  <j>  dll’,  (22) 

where  r/(x,  t)  is  the  surface  elevation.  For  the  energy  computations,  the  complex  velocity  potential 
is  estimated  by  integrating  the  complex  velocity  along  the  free  surface.  We  use  80  elements  and 


9 


a  time  step  of  0.2.  The  time  step  is  sufficiently  small  such  that  time  discretization  errors  are 
negligible  compared  to  spatial  discretization  errors. 

The  potential  method  is  more  accurate  for  short  times  than  the  velocity  method  for  total 
energy  conservation,  while  both  methods  show  good  agreement  for  the  other  physical  and  error 
quantities.  After  two  wave  periods,  however,  the  initial  energy  level  is  recovered  in  the  velocity 
formulation.  In  the  velocity  formulation,  the  complex  velocity  must  be  integrated  numerically  in 
order  to  determine  the  kinetic  energy.  VVe  can  integrate  k  analytically  in  a  panel,  because  linear 
variation  of  k  in  the  mapped  domain  is  assumed  in  the  panel.  Other  integrating  schemes  using 
B-spline  interpolation  resulted  in  a  maximum  energy  loss  which  is  twice  as  large.  Hence,  we 
expect  that  high-order  discretization  of  the  Cauchy’s  integral  will  improve  energy  conservation. 

The  present  method  is  very  sensitive  to  the  computation  of  the  tangent  angle.  As  the  non¬ 
linear  free  surface  wave  evolves,  the  local  curvature  can  be  very  large,  so  the  reliability  of  the 
present  scheme  can  deteriorate.  This  difficulty  may  be  resolved  bv  adopting  a  carefully  chosen 
remeshing  scheme. 

Figure  4  shows  the  results  of  a  convergence  test  for  Em  and  the  energy  conservation  criteria, 
Ee,  at  fixed  time  t=  1 .0.  This  test  increases  the  number  of  elements  N  from  20  to  80  while  holding 
the  time  step  fixed  at  At  =  0.1.  Ee  is  defined  as 

Et  ={E{t)  -E(0))/E(0).  (23) 

Both  methods  show  very  good  (and  identical)  convergence  for  the  mass  conservation  criteria, 
but  poorer  convergence  for  energy  conservation.  These  examples  also  show  that  the  velocity 
formulation  requires  slightly  more  computational  time  (2  percent)  than  the  potential  method, 
due  to  the  need  to  calculate  the  tangent  angle  and  the  material  derivative  of  the  tangent  angle. 


10 


4.2  Vortex  Pair  and  Wall  Interaction 

As  another  example,  the  convection  of  a  vortex  pair  approaching  a  solid  wall  is  examined  to 
test  the  accuracy  of  the  numerical  algorithm.  The  exact  solutions  for  the  induced  velocity  and 
the  path  of  the  vortex  pair  are  given  in  Milne-Thomson  [13].  Since  the  tangent  angles  on  the 
boundary  contour  (solid  wall)  are  fixed  in  this  example,  the  errors  in  the  approximation  of  the 
tangent  angle  (as  in  section  4.1)  are  eliminated. 

The  computational  domain  is  reduced  to  a  finite  one  by  solving  a  periodic  problem  with  a 
period  sufficiently  long  so  that  the  periodic  boundary  effects  are  negligible.  The  conformal  map¬ 
ping  of  the  first  example  is  again  used.  The  initial  location  and  strength  of  the  vortex  pair  are 
given  in  Table  2  along  with  the  imposed  period. 

Figure  5  shows  the  initial  error  in  tangential  velocity  using  40  equally  spaced  elements.  The 
error  is  normalized  by  the  solution  maximum.  It  shows  that  the  present  scheme  gives  an  average 
error  norm  of  one  part  in  10~6. 

The  convergence  rate  of  the  tangential  velocity  is  examined  by  plotting  the  root  mean  square 
error  (E2)  as  a  function  of  the  number  of  elements  (N)  using  log-log  scales  shown  in  Figure  6. 
Surprisingly,  the  present  scheme  is  fourth-order  accurate  in  this  example.  Usually,  the  results 
of  the  linear  approximation  are  second-order  accurate.  This  higher  convergence  rate  can  be  at¬ 
tributed  to  the  use  of  an  overdetermined  system  on  a  circular  (after  conformal  mapping)  domain 

[9]- 

\ 

l 

On  the  other  hand,  the  tangential  velocity  by  the  potential  formulation  is  only  second-order 
accurate.  This  problem  can  also  be  solved  by  the  potential  formulation,  because  the  influence  of 

! 

the  vortices  on  the  wall  boundary  condition  (rj>)  is  uniquely  given.  Even  if  the  velocity  potential 
1  calculated  by  the  potential  formulation  were  fourth-order  accurate,  the  accuracy  would  be  low- 
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ered  by  numerical  differentiation  (second-order  central  difference)  to  get  the  disturbance  velocity 
components.  The  total  velocity  field  in  this  case  is  calculated  as  the  sum  of  the  numerically 
calculated  disturbance  velocity  and  of  the  direct  contribution  from  the  field  vortices. 

The  computed  path  of  the  vortex  pair  in  time  is  given  in  Figure  7  along  with  the  analytic 
solution.  This  shows  that  the  time  stepping  method  used  in  the  present  study  is  sufficiently 
accurate  to  follow  the  convected  motion  of  the  vortex  pair. 

4.3  Vortex  Pair  and  Free  Surface  Interaction 

Our  final  case  examines  a  vortex  pair  interacting  with  a  free  surface.  Initially  the  vortex  pair 
is  introduced  into  the  fluid  with  the  free  surface  undisturbed.  Except  for  the  vortex  pair,  the 
fluid  is  assumed  to  be  irrotational. 

In  a  potential  formulation,  the  free  surface  boundary  condition  is  given  by  the  velocity  po¬ 
tential.  But  the  velocity  potential  due  to  vortices  is  multivalued,  as  we  discussed  in  section  2. 
To  treat  this  multivalued  potential  in  the  numerical  computation,  a  carefully  determined  logical 
statement  must  be  used.  For  an  arbitrary  boundary  contour,  it  is  difficult  to  devise  an  efficient 
logical  statement  for  uniqueness.  Moreover,  when  the  vectorization  of  the  computer  code  is  de¬ 
sirable,  we  cannot  use  a  logical  statement.  Hence,  only  the  result  of  a  velocity  formulation  is 
shown  here. 

In  order  to  reduce  the  computational  domain  to  a  finite  one,  damping  layers  are  introduced 
on  both  ends  of  computational  domain  following  Telste’s  example  [14].  The  initial  conditions  and 
strength  of  the  vortex  pair  are  given  in  Table  3.  Initially,  nodes  are  distributed  using  a  cosine 
spacing  so  that  more  material  points  are  located  near  the  origin.  This  ensures  adequate  spatial 
resolution  near  the  origin  for  a  short  time,  since  the  continuous  upwelling  of  the  free  surface 
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convects  nodes  away  from  the  counter-rotating  vortices. 

In  Figure  8,  the  free  surface  elevations  at  t  =  13.4164  are  given.  The  results  of  Telste  [14], 
using  Bilker’s  [3]  generalized  vortex  method  with  damping  layers,  are  in  good  agreement  with 
our’s,  although  we  use  a  much  smaller  number  of  nodes  (N  =  81  compared  to  641).  Unlike 
the  example  in  section  4.1,  the  Lagrangian  movement  of  the  nodes  grows  rapidly  as  the  vortex 
pair  approachs  the  free  surface.  Hence,  a  remeshing  scheme  of  free  surface  nodes  is  needed  for  a 
reliable  solutions  at  longer  times. 

5  Conclusions 

A  velocity  form  of  the  complex  variable  boundary  element  method  was  developed  for  po¬ 
tential  flow  by  applying  Cauchy’s  integral  theorem  to  the  complex  velocity.  A  new  form  of  the 
full  nonlinear  dynamic  free  surface  boundary  condition  was  derived  as  a  function  of  the  normal 
and  tangential  velocity  components  and  successfully  applied  to  three  examples. 

The  main  advantage  of  the  present  method  over  the  conventional  potential  method  is  that 
it  can  more  easily  solve  for  flow  with  field  (free)  singularities  such  as  vortices  and  sources.  The 
present  method  also  gives  the  velocity  field  directly  without  the  need  for  numerical  differentiation. 
However,  it  is  very  sensitive  to  the  definition  of  the  local  tangent  angle  of  the  free  surface  so  that 
remeshing  and  filtering  are  necessary  to  assure  a  reliable,  stable  solution.  Specifically,  the  velocity 
formulation  does  not  conserve  energy  as  accurately  as  the  potential  formulation  for  free-surface 
problems  without  vorticity.  This  could  be  due  to  errors  in  determining  the  local  tangent  angle  or 
in  determining  the  velocity  potential  by  numerical  integration,  which  is  required  for  the  kinetic 
energy  calculation. 


Acknowledgements 


This  work  was  partially  supported  by  Naval  Research  Laboratory  Contract  No.  N00014-85-K- 
Z019,  ONR  Contract  N000 14-87-0509,  and  the  Program  in  Ship  Hydrodynamics  at  The  Univer¬ 
sity  of  Michigan,  funded  by  the  University  Research  Initiative  of  the  Office  of  Naval  Research, 
Contract  No.  N000 184-86- K-0684.  The  authors  acknowledge  G.  Trygvasson  for  his  suggestions. 


14 


References 


1.  M.S.  Longuet-Higgins  and  C.D.  Cokelet,  “The  deformation  of  steep  surface  waves  on  water: 
I.  A  Numerical  Method  of  Computation,”  Proc.  R.  Soc.  London,  AS50,  1-26  (1976). 

2.  T.  Vinje  and  P.  Brevig,  “Numerical  calculations  of  forces  from  breaking  waves,”  Hydrody¬ 
namics  in  Ocean  Engineering,  Norwegian  Inst.  Tech.,  547-565  (1981). 

3.  G.R.  Baker,  D.I.  Meiron  and  S.A.  Orszag,  “Generalized  vortex  methods  for  free-surface 
flow  problems,”  J.  Fluid  Meek.  1SS,  477-501  (1982). 

4.  J.W.  Dold  and  D.H.  Peregrine,  “Steep  unsteady  water  waves:  An  efficient  computational 
scheme,”  19th  Int.  Conf.  on  Coastal  Engineering,  Houston  (1984). 

5.  E.Martensen,  ARMA,  3,  235-270  (1959)  (in  German). 

6.  R.I.  Lewis  and  P.G.  Ryan,  J.  Meek.  Eng.  Sci.,  14,  280  (1972). 

7.  A.  Leonard,  “Vortex  methods  for  flow  simulation,”  J.  Comp.  Phys.  37,  289-335  (1980). 

8.  G.R.  Baker,  D.I.  Meiron,  and  S.A.  Orszag,  “  Generalized  Vortex  Methods  for  Free-Surface 
Flow  Problems;  II  Radiating  Waves,”  submitted  for  publication  (1986). 

9.  W.W.  Schultz  and  S.W.  Hong,  “Solution  of  potential  problems  using  an  overdetermined 
complex  boundary  integral  method,”  submitted  for  publication.  (1988) 

10.  S.W.  Hong,  G.  Tryggvason  and  W.W.  Schultz,  “A  note  on  the  formulation  of  free  surface 
problems.  An  integral  equation  for  velocity,”  submitted  for  publication.  (1987) 

11.  R.  Fletcher  and  C.M.  Reeves,  Computer  Journal,  7,  149-154  (1964). 

12.  W.W.  Schultz,  Proceedings  of  11th  IMACS  World  Congress  S,  Oslo,  Norway  (1985). 

13.  L.M.  Milne  and  C.B.E.  Thomson  ,  Theoretical  Hydrodynamics,  The  University  Press,  Glas¬ 
gow.  (1968) 


15 


14.  J.G.  Telste,  “Potential  flow  about  the  counter  rotating  vortices  approaching  a  free  surface,” 
submitted  for  publication.  (1987) 


16 


Table  1:  Initial  Conditions  and  Other  Parameters 
for  Nonlinear  Free  Surface  Problem 


Imposed  Period 

2n 

Domain  of  Problem 

0  <  x  <  2tt 

Length  Scale  L 

Period/27r 

Velocity  Scale 

(gL)* 

Surface  Elevation 

y(x)  =  asinx 

Complex  Potential 

/?  =  a  e~* 

Tangential  Velocity 

a  e5'(sinx  —  acos2x)/(l  +  a2cos2x)i 

Amplitude 

o 

II 
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Table  2:  Initial  Conditions  and  Other  Parameters 


9 

I 

I 


for  Wall- Vortex  Problem 


Imposed  Period 

2ir 

Domain  of  problem 

—  TT  <  X  <  7 r 

Length  Scale  L 

Period/27r 

Velocity  Scale 

Time  Scale 

miwm 

Location  of  Vortex 

Strength  of  Vortex 

(100),  (-100) 
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Table  3:  Initial  Conditions  and  Other  Parameters 
for  Free  Surface- Vortex  Problem 


Domain  of  problem 

-67  <  x  <  67 

Damping  Layer 

47  <  |xj  <67 

Length  sacle  L 

Initial  Distance  of  Two  Vortices 

Velocity  Scale 

Time  Scale 

{L/g)* 

Location  of  Vortex 

(0.5, -5.0), (-0.5, -5.0) 

Strength  of  Vortex 

(2.236068), (-2.23606S) 

LIST  OF  FIGURES 


Figure  1:  Domain  and  Bounding  Surfaces 

Figure  2:  Comparison  of  Errors  in  Velocity  Between  Two  Methods 
(N=40,  t=0.0) 

Figure  3:  Variation  of  Ymax,  Fm,n,  Em  and  E(t )  in  Time 
(N=80,  At  =  0.2) 

Figure  4:  Results  of  Convergence  Test  for  Em,  and  Ee 
(At  =  0.1,  t— 1.0) 

Figure  5:  Error  in  Tangential  Velocity  for  Wall- Vortex  Problem 
(N=40,t=0.0) 

Figure  6:  E2  Errors  in  Tangential  Velocity  for  Wall- Vortex  Problem  (t=0.0) 

Figure  7:  Locus  of  Vortex  Motion 

Figure  8:  Free  Surface  Elevation  at  t= 13.4164 


Figure  1:  Domain  and  Bounding  Surfaces 


U  ERR 


(b)  y  velocity 


Fig“re  ^  (N0“oV  ErrOTS  Vel0d*  Two  Methods 


V  .  vvv  ^ 

0  . 


20 .  40 .  60 .  80 . 

N 


Figure  4:  Results  of  Convergence  Test  for  Em,  and  E, 
(At  =  0.1,  t=1.0) 


0.000005 


0 . 000000 


-0 . 000005 
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Figure  6:  E?  Errors  in  Tangential  Velocity  for  Wall- Vortex  Problem  (t  0.0) 
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